#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>
#include <fstream>
#define MAXBIT 30
#define SCRAM 3
#define INIT 1 

namespace sobol{

extern  "C"{
void sobol_(double *x, int *n, int *d, double *last, int *LL, int *count, int *SV, int *scrambling, int *seed, int *init, int *trans);
}	
//SV of size dimen*MAXBIT=30
	//last:1*d
	//x=n*d
	//C     scrambling     - INITIALIZATION FLAG
	//C     0 - NO SCRAMBLING
	//C     1 - OWEN TYPE SCRAMBLING
	//C     2 - FAURE-TEZUKA TYPE SCRAMBLING
	//C     3 - OWEN + FAURE-TEZUKA TYPE SCRAMBLING
	// TRANS - FLAG, 0 FOR UNIFORM, 1 FOR NORMAL DISTRIBUTION
/* Allocates a double matrix with subscript range m[0...nr-1][0...nc-1] */

int getSeed()
{
    std::ifstream rand("/dev/urandom");
    char tmp[sizeof(int)];
    rand.read(tmp,sizeof(int));
    rand.close();
    int* number = reinterpret_cast<int*>(tmp);
    return (*number);
}

double ** matrix(int nr, int nc)
{
	 int ii;
	 double **mat;
	 /* allocate pointers to rows */
	 mat=(double**)calloc(nr,sizeof(double*));
	 /* allocate rows and set pointers to them */
	 if (!mat) printf("%s","Allocation failure 1 in matrix()");
	 for (ii=0; ii<nr; ii++) {
	      mat[ii]= (double*)calloc(nc,sizeof(double));
	      if (!mat[ii]) printf("%s","Allocation failure 2 in matrix()");
	}
	return(mat);
}

template<typename T>
T* vector(int n)
{
	  T *v;
	  v=(T*)malloc(n*sizeof(T));
	  if (!v) printf("%s","Allocation failure in vector()");
	  return(v);
}
template<typename T>
void free_vector(T *v)
{
  free(v);
}

/* Free a matrix */
void free_matrix(double **mat, int nr)
{
  int ii;
  for (ii=(nr-1);ii>=0;ii--) free((char*)(mat[ii]));
  free((char*)(mat));
}
class Owenscramble{

	public:
		Owenscramble(void){
			boost::random::uniform_01<> U;
			runif=new RandomG::Random<boost::random::uniform_01<> >(U);
		}
double*  operator()(int n, int d, int Gaussian,double *x)
{
	//double *x=vector<double>(d*n);
	double *last=new double[d];//vector<double>(d);
	int *SV=new int[MAXBIT*n];//vector<int>(MAXBIT*n);
	int scram=SCRAM;
	int seed=getSeed();
	int init=INIT;
	int ll=1;
	int count=1;
	int k=0,i=0,j=0;
	for(i=0;i<d;i++)
	{
		last[i]=(*runif)();
		
	}
	for(j=0;j<MAXBIT*n;j++)
	{
		SV[j]=j;
	}
	
	sobol_(x,&n,&d,last,&ll,&count,SV,&scram,&seed,&init,&Gaussian);
	delete[] last;
	delete[] SV;
	//free_vector<double>(x);
	//free_vector<double>(last);
	//free_vector<int>(SV);
//dont forget to free
	return x;
}
private:
	double *current;
	RandomG::Random<boost::random::uniform_01<> > *runif;
};
}//end sobol namepace
